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We present a detailed analysis of the mechanism of the primary charge separation process in bacterial photosynthesis using 
real-time path integrals. Direct computer simulations as well as an approximate analytical theory have been employed to map 
out the dynamics of the charge separation process in many regions of the parameter space relevant to bacterial photosynthesis. 
Two distinct parameter regions, one characteristic of sequential transfer and the other characteristic of superexchange, have 
been found to yield charge separation dynamics in agreement with experiments. Nonadiabatic theory provides accurate rate 
estimates for low-lying and very high-lying bacteriochlorophyll state energies, but it breaks down in between these two regimes. 
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I. INTRODUCTION 

The problem of bacterial photosynthesis has at- 
tracted a lot of recent interest since the structures 
of the photosynthetic reaction center (RC) in purple 
bacteria Rhodopseudomonas viridis and Rhodobacterias 
sphaeroides have been determined jl] . Much research ef- 
fort is now focused on understanding the relationship be- 
tween the function of the RC and its structure. One fun- 
damental theoretical question concerns the actual mech- 
anism of the primary electron transfer (ET) process in 
the RC, and two possible mechanisms (or a combina- 
tion of both) have emerged out of the recent theoretical 
work 1^-^. The first is an incoherent two-step mecha- 
nism where the charge separation involves a sequential 
transfer from the excited special pair (P*) via an inter- 
mediate bacteriochlorophyll monomer (B) to the bacte- 
riopheophytin (H). The other is a coherent one-step su- 
perexchange mechanism, with P'*'B~ acting only as a vir- 
tual intermediate. 

From experimental studies to date, one cannot eas- 
ily rule out either possibility. In fact, conflicting inter- 
pretations have been derived from different experiments. 
The detection of transient population in the P+B^ state 
has been considered by many to be the key experimental 
evidence that can differentiate between the two mecha- 
nisms. Many transient absorption spectroscopic experi- 
ments [§,0 could not detect bleaching in the B~ band, 
lending support to the coherent superexchange mecha- 
nism. On the other hand, some other experiments [ pl]JT2| 
point toward a two-step process since a fast second rate 
constant corresponding to the P+B~H P+BH~ transi- 
tion was detected. The picture emerging from a compar- 
ison of the molecular dynamics (MD) simulations for the 
RC carried out so far p^-p^ is equally murky. Different 
groups arrive at conflicting conclusions. 



Despite the emphasis placed on the experimental de- 
tection of population in the P"'"B~ band, we have shown 
in a recent real-time path integral simulation [ p^ that ob- 
servation of transient population in P+B^ cannot deflni- 
tively rule out either mechanism. In fact, significant pop- 
ulation in the P+B" state (up to 20%) can be observed 
in both the sequential and the superexchange regimes in 
the simulations. This suggests that transient population 
on the P+B^ state, even if detectable, may not really 
be useful for differentiating between the two mechanisms 
0. In fact, almost all experimentally observed dynami- 
cal features of the RC can be reproduced nicely in either 
regime. 

In this Letter, we report a detailed path integral anal- 
ysis of the mechanism of the primary charge separation 
process. We have computed the dynamics of the primary 
ET (i.e., the time-dependence of the occupation probabil- 
ities on the three relevant chromophorcs) as a function of 
the P+B^ energy for many combinations of the electronic 
coupling and the reorganization energy. The primary ob- 
jective is to determine the parameter values that are re- 
quired for the dynamics to be consistent with character- 
istic experimental observations, and hence to provide a 
map of the various possible parameter regions in which 
the RC could operate. Our previous simulations have 
been carried out for a relatively small energy range for the 
P+B~ state (between —666 cm~^ and -1-666 cm~^ rela- 
tive to the P* state). In this Letter, we present new re- 
sults for a much larger range of energies, up to 3000 cm~^ 
above and 1333 cm~^ below the P* state. This spans the 
region usually considered to be relevant for the RC, ex- 
cept for the region with extremely high P"'"B~ en ergy (a 
value of -1-8000 cm^^ has been suggested in Ref. |13|] ). 
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II. THREE-STATE SPIN-BOSON MODEL FOR 
THE RC 

The model for our studies has been discussed previ- 
ously in great detail [^,^ 16 1. For the sake of complete- 
ness, we will summarize some of the relevant features 
here. For a model of the RC, we consider Hamiltonians 
of the form 



H = Ha + HE 



(1) 



where Hq describes the bare three-state system, Hb 
corresponds to a Gaussian bath describing the protein 
environment, and Hj is a bilinear system-environment 
coupling describing the interaction of the electric dipole 
moment with the bath polarization. The states 1, 2, 
and 3 correspond to the electronic configurations P*BH, 
P+B"!! and P+BH~, respectively, and we allow for ar- 
bitrary binding energies of these states. The bare three- 
state Hamiltonian is 



El -Ku -Ki3 
Ho = I -Ki2 E2 -K23 
-K13 —K23 E3 



(2) 



where the electronic coupling (tunnel matrix element) 
between states i and j is denoted by Kij . 

This spin-boson model is an idealized model for the 
RC. Its limitations have been discussed in Ref. |jl^. The 
most significant aspect of this model, however, is its abil- 
ity to capture the effects of the low-frequency protein 
modes on the tunneling electron from first principles p7| . 
In this sense, the protein casing in the RC presents it- 
self to the electron as a Gaussian fluctuating polariza- 
tion field with a continuous spectral density. Previous 
MD simulations have shown that this spectral density 
is smooth, featureless, and of Ohmic form with a char- 
acteristic frequency lUc of about 160 cm~^. They have 
also confirmed the assumption of Gaussian statistics for 
the bath polarization since the free energy curves were 
always found to be strictly parabolic. In the classical 
(high temperature) limit, the spectral density has a sim- 
ple relationship with the familiar reorganization energy 
Ai3 between states 1 and 3, and one can characterize the 
Ohmic spectral density completely by the two parame- 
ters ujc and A13. We fix a number of parameters in our 
model according to reasonable experimental and theoret- 
ical estimates |l6| : 

1. With the energy scale fixed by setting Ei = 0, the 
energy of the P+H" state E3 is set at —2000 cm~^. 

2. The bath frequency Wc is set at 166 cm^^ according 
to MD results. 

3. The ratio K23/K12 is set at 4, and K13 is assumed 
to be negligible. All electronic couplings are as- 
sumed to be independent of the actual protein con- 
figuration (Condon approximation). 



The technical details of our real-time quantum Monte 
Carlo (QMC) simulations have been reported in great 
detail elsewhere and will not be repeated here. 
This numerical technique allows for a computation of 
the time-dependent occupation probabilities on the elec- 
tronic states up to a certain time limit. QMC simula- 
tions cannot be carried out to very long times due to the 
well-known dynamical sign problem. However, the simu- 
lations reported here have been carried out up to a few 
picoseconds, the same timescale as the experimentally 
determined ET rate. In these simulations, we search for 
parameter regions in which the charge separation pro- 
ceeds in qualitative agreement with experimental obser- 
vations. We consider the following to be key experimental 
characteristics: (1) the P+B~ population (P2) should be 
small throughout ( ^ 20%); (2) the charge separation rate 
should be about 3 ps at room temperature; and (3) the 
rate should increase about twofold at cryogenic temper- 
atures. For each value of the P+B^ state energy E2, we 
attempt to find a combination of K12 and A13 that would 
yield dynamics with these experimental characteristics. 



III. QMC SIMULATION RESULTS AND PHASE 
DIAGRAM 

From the QMC simulations, we have compiled a "phase 
diagram" for the RC, which is a map of the various re- 
gions in parameter space which generate dynamics with 
all the key experimental characteristics. This phase dia- 
gram is shown in Fig. 1, with the proper values of K12 
correlated with the energy of the P+B~ state i?2- On the 
phase diagram, we have also indicated the values of K12 
predicted by conventional nonadiabatic theory to give a 
room-temperature ET rate of 3 ps (without demanding 
that they also conform to the inverse temperature depen- 
dence) . 

From the phase diagram, there are clearly two distinct 
regions which yield charge separation dynamics with the 
correct experimental characteristics. The first is char- 
acterized by a P+B~ state energy E2 lower than about 
—600 cm^^ and down to about —1300 cm^^. Within 
this region, the dynamics largely agrees with predic- 
tions from conventional nonadiabatic theory for the se- 
quential mechanism. Simulation results for the transient 
populations on the three electronic states are shown in 
Fig. 2 at two temperatures for E2 = —666 cm~^ and 
— 1333 cm~^. The 1 — > 3 transfer is efficient in both 
cases, with P2 remaining below 15% throughout. For 
E2 = —666 cm^^, the charge separation rate increases 
from 3.0 ps at 298 K to 2.0 ps at 140 K. This should be 
contrasted with E2 = —1333 cm"-*^, for which the rate 
increases from 3.3 ps at 298 K to only 2.7 ps at 140 K. 
Although the ET rate still increases with lower tempera- 
ture for E2 — —1333 cm~^, this temperature-dependence 
is somewhat weaker than that observed in experiments. 
For E2 much below —1333 cm^^, the temperature depen- 
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dence is not strong enough to be consistent with exper- 
iments. Therefore, in the phase diagram, the sequential 
region terminates around E2 — —1333 cm^^. 
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FIG. 1. Phase diagram of the RC compiled from QMC 
simulation data, correlating 7^12 with E2 (thick lines). The 
proper reorganization energies are indicated on the phase dia- 
gram for those parameter sets that yield dynamics consistent 
with experiments (open circles). Nonadiabatic predictions for 
the 7^12 yielding a room-temperature rate of 3 ps are also 
shown for comparison (dashed curve). 

The second region is characterized by a P+B^ state en- 
ergy E2 higher than -1-666 cm~^ (we have performed sim- 
ulations for E2 up to -1-3000 cm~^). For E2 = +666 cm~^ 
and -1-1333 cm~^, conventional nonadiabatic superex- 
change theory does not give reliable predictions for the 
dynamics. The transfer mechanism in this region is pre- 
dominantly superexchange, although the precise rates are 
not well described by nonadiabatic theory. As the energy 
of the P+B~ state moves up to E2 = -1-2000 cm~^ and be- 
yond, the dynamics becomes increasingly nonadiabatic. 
We mention in passing that in this limit the dynamical 
sign problem is quite severe, and one needs to sample 
extremely long Monte Carlo trajectories rendering QMC 
simulations for E2 much higher than -1-3000 cm~^ very 
costly. Nevertheless, our present data show that for suf- 
ficiently high-lying P+B^ state, the conventional golden 
rule superexchange formula becomes increasingly accu- 
rate, making QMC simulations less valuable for E2 > 
+3000 cm"^. For E2 < -1-2000 cm"\ however, nonadi- 
abatic theory seems to severely underestimate the mag- 



nitude of the electronic couplings required to achieve su- 
perexchange rates consistent with experiments. 
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FIG. 2. QMC data for the ET dynamics with low-lying 
P+B" state and K12 = 28 cm"^ for E2 = -666 cm"^ 
(Ai3= 4000 cm^i) and E2 = -1333 cm'^ (Ai3= 5300 cm'^). 
The dashed (dotted, solid) curve gives the occupation proba- 
bility on state 1 (2, 3). 



Simulation results for the transient populations on 
the three electronic states are shown in Fig. 3 at two 
temperatures for E2 = +666 cm~^, +1333 cm~^ and 
-f2000 cm~^. For all three energies, the charge separa- 
tion rates at 298 K are about 3 ps. These rates increase 
to about 2.2 to 2.5 ps at 140 K. As was pointed out in 
our previous work | |l6[ |, the dynamics in this regime is 
not simply monoexponential. There is a fast short-time 
component for the first w 0.5 ps, followed by a slower 
component. The fast component always obeys the ex- 
perimentally observed inverse temperature dependence, 
but the rate and the proportion of the slower component 
vary nonmonotonically with temperature, as measured in 
many mutant RCs p8| . The physical reason for this phe- 
nomenon has been given by Gehlen et al. |p^ , and our 
simulations provide direct numerical support for their ar- 
guments. 

The phase diagram in Fig. 1 also reveals another in- 
teresting aspect of the superexchange regime: the RC 
seems to be able to operate under a wide range of P+B~ 
state energies with the same conditions, i.e., with ap- 
proximately the same electronic coupling K12 (about 140 
to 150 cm~^) and the same reorganization energy A13 
(about 2000 cm^^). This means that within the superex- 
change regime, the system is robust against variation in 
the P+B" state energy. The sequential region is not 
as robust. Although approximately the same value of 
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Ki2 (about 28 cm"-'^) would work for E2 between about 
— 1333 cm~^ and —666 cni~^, the optimal value of the 
reorganization energy A13 changes from 4000 cm~^ at E2 
= -666 cm-i to 5300 cir-^ at £'2 = -1333 cm^^. 
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FIG. 3. QMC data for high-lying P+B~ state with 
Ai3= 2000 cm~^. The couplings are Kn = 140 cm"^ 
for El = +666, +1333 cm-\ and K12 = 150 cm"^ for 
E2 = +2000 cm-^ 

Finally, it is important to stress that the RC 
does not seem to operate via a combined sequen- 
tial/superexchange mechanism as proposed in Refs. 
|^,|o|. Our simulations indicate clearly that it is cither 
a pure sequential or a pure superexchange transfer. This 
is indicated by the grey region in Fig. 1, which we refer 
to as the "gap" region. We emphasize that this does not 
imply that in this region it is impossible to achieve a fast 
ET rate. On the contrary, it is possible to achieve a 3 ps 
transfer at room temperature with only minor population 
accumulation on the P+B~ state. However, we have not 
found any parameter set in this region that can simul- 
taneously satisfy all the experimentally observed charac- 
teristics both at high and low temperatures pq] . 



IV. NONINTERACTING-CLUSTER 
APPROXIMATION 

By invoking a simple and reasonable approximation 
for the formally exact path-integral expressions, we 
have previously derived a set of nonlocal master equa- 
tions for the time-dependent occupation probabilities . 
Our "noninteracting-cluster approximation" (NICA) ex- 
presses the time-dependent transition rates between the 
various sites in terms of cluster functions Fy . These clus- 
ter functions give the amplitude of all paths going from 
diagonal state i to diagonal state j of the reduced den- 
sity matrix without touching the diagonal in between. 
Since each hop on the 3x3 lattice representing the pos- 
sible states of the reduced density matrix gives a factor 
HKij, the rates can be expressed as a power series in 
the electronic couplings. 

By confining this expansion to the lowest-order terms 
and taking the long-time limit (which makes the master 
equations local in time), one obtains the nonadiabatic 
rate expressions referred to in Sec. III. The rates and 
F23 (plus backward rates) describing a sequential transfer 
are just golden rule rates. Furthermore, the lowest-order 
cluster ~ ^12^23 gives the nonadiabatic superex- 
change rate. There is only one pathway contributing to 
this fourth-order cluster (plus the complex conjugate), 
namely the one going along the edge of the 3x3 lattice. 
The resulting rate formula can be further simplified in the 
classical limit by taking a short-time approximation for 
the bath kernel, or for a high-lying state where one can 
derive the conventional golden rule rate for the emerging 
two-state system spanned by states 1 and 3 These 
golden rule rates are consistent with the QMC results for 
very high P+B~ state energy (£'2 si +2000 cm~^) as well 
as for low-lying P+B" states {E2 ^ — 666 cm^^). But in 
the intermediate region +666 cm~^ ^ £2 ^ +2000 cm^^, 
nonadiabatic theory breaks down (see Fig. 1). 

We will now study the origin of the apparent dis- 
crepancy between nonadiabatic theory and the QMC re- 
sults in this intermediate region. To that purpose, we 
have computed the next-order corrections to the superex- 
change cluster, r^g"*. Here one has to consider the two 
pathways with six hops which go from diagonal state 1 
to state 3 without hitting the diagonal. In the classical 
limit, we obtain for the activationless situation consid- 
ered in the simulations (£3 = — A13) 



(6) 
13 
(4) 
13 
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2(3F* 



l + exp{(3F*) 



(3) 



where 5E = E2 + A13/4 and F* = {SE)'^/\i3. For 
a high-lying intermediate state, the bracket gives unity 
because f3F* ^ 1. Therefore the sixth-order correc- 
tion decreases the nonadiabatic estimate by ~ 20% for 
E2 = +666 cm~^. With the short-time bath kernel, it 
would actually be possible to sum the whole power series 
and thus to obtain the full classical NICA superexchange 
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rate without stopping at sixth order. However, we have 
not done so because this short-time approximation is in- 
appropriate for high-order diagrams. In any ease, it is 
expected that these higher-order contributions lead to a 
reduction of the nonadiabatic superexchange rate. More- 
over, we also expect that adiabatic corrections beyond 
NICA can lead to substantial renormalizations for large 
electronic couplings. 

To conclude, we have given a path-integral analysis 
of the primary electron transfer step in bacterial pho- 
tosynthesis. The mechanism appears to be either a se- 
quential or a superexchange transfer but not a combined 
one. Whereas conventional nonadiabatic theory appears 
to be accurate for both the sequential and the ultimate 
superexchange regime, adiabatic corrections reduce the 
superexchange rate substantially for intermediate ener- 
gies of the P+B^ state. 
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